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Abstract 

An implicit, finite-difference procedure for 
numerically solving viscous incompressible flows is 
presented. The pressure-field solution is based on 
the pseudocompressibility method in which a time- 
derivative pressure term is introduced into the mass- 
conservation equation to form a set of hyperbolic equa- 
tions. The pressure-wave propagation and the spread- 
ing of the viscous effect is investigated using simple test 
problems. Computed results for external and internal 
flows are presented to verify the present method which 
has proved to be very robust in simulating incompres- 
sible flows. 

I. Introduction 

Incompressible flow phenomena are frequently 
encountered in many engineering applications, espe- 
cially, in hydrodynamics and in certain classes of 
aerodynamic problems such as dynamic stall and low- 
speed wind-tunnel test problems. For most two- 
dimensional flow simulations, computer time and 
memory requirements are not major limiting factors, 
and various numerical techniques have been imple- 
mented quite successfully. For example, a stream 
function-voirticity formulation is frequently used for 
solving two-dimensional, viscous, incompressible flow 
problems (for example, see ref. 1). This approach 
conveniently eliminates the pressure terms. The three- 
dimensional extension of this method, however, is 
not straightforward. Various other three-dimensional 
viscous-flow solvers have been developed, mainly, 
for compressible flow (for example, see refs. 2-5). 
Implementing these procedures for simulating viscous 
incompressible flows is not efficient and is generally not 
recommended. Therefore, it is desirable to develop an 
efficient computational method for solving viscous in- 
compressible flows. For convenience in applying the 
method to three-dimensional problems, primitive vari- 
ables, namely, the pressure and velocities, are used in 
the present work. 

One of the major difficulties in solving incom- 
pressible flows that use primitive variables is caused 
by the presence of the pressure term which is used 
as a mapping parameter to obtain a divergence-free 
velocity field. One method of solving for pressure is 
the Poissons equation approach developed by Harlow 


and Welch [6], which has been used frequently, mostly 
in conjunction with explicit methods. The usual com- 
putational procedure is to choose the pressure field 
such that continuity is satisfied at the next time- 
level, so that the new flow field will be divergence-free. 
This procedure normally requires a relaxation scheme 
iterating on pressure until the divergence-free condi- 
tion is reasonably satisfied. Because this approach 
can be very time consuming, the computing time re- 
quired for simulating three-dimensional flows has been 
prohibitively large. To accelerate the pressure-field 
solution and alleviate the drawback associated with 
the Poissons equation approach, Chorin [7] proposed 
the use of artificial compressibility in solving the con- 
tinuity equation. A similar method was adopted by 
Steger and Kutler [8] using an implicit approximate- 
factorization scheme due to Beam and Warming [9]. 
To implement the implicit time-differencing, a hyper- 
bolic time-dependent system of equations is fabricated 
by adding a time-derivative of the pressure term to the 
mass-conservation equation. These hyperbolic equa- 
tions possess characteristics that are not present in 
the usual Poisson’s equation for pressure. This ap- 
proach has been applied to simulate laminar, incom- 
pressible flow within liquid- filled shells [10], Based 
on this procedure, a pseudocompressible method has 
been developed for solving three-dimensional, viscous, 
incompressible flow problems cast in generalized cur- 
vilinear coordinates [11,12]. The purpose of the present 
paper is to show salient features of the pseudocompres- 
sible approach using numerical experiments. 

The method is described in section n, the numeri- 
cal algorithm is explained in section HI, and computa- 
tional experiments are presented in section IV to show 
the physical characteristics of the present method. 

II. Description of the Method 
Incompressible Navier-Stokes Equations 

Unsteady, three-dimensional, viscous, incom- 
pressible flow with constant density is governed by the 
following Navier-Stokes equations, written in tensor 
notation: 

S = ° (Xa) 

duj dujUj _ dp dm . . 

dt i ' dxj dxi ' T dxj { ' 

Here, t is time; Xi are the Cartesian coordinates; u t are 
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corresponding velocity components; p is the pressure; 
and Tij is the viscous stress tensor. The viscous stress 
tensor can be written in the following form: 

T\j 2 V S ij Rij (lc) 


The present method is primarily designed for ob- 
taining steady-state solutions efficiently, and our inter- 
est here is limited to the development of the steady- 
state algorithm. 


The strain rate tensor is defined by 


Sa 


1 ,duj 

2 ( dXj 


+ 


duj 

dxi 


(Id) 


Here, ‘is the Reynolds stresses, and u is the 

coefficient of viscosity. Since the expression for 
the Reynolds-stress terms depends on the turbulence 
closure model used, the accuracy of the turbulent flow 
problems will naturally depend on the level of the 
model and on experimental input used for calibrat- 
ing the model. Because economy is still an overriding 
factor in many engineering applications, turbulence is 
often simulated by an algebraic eddy-viscosity model, 
using a constitutive equation involving a " mixing 
length” as a measure of the turbulence length scale. 
In the present study, however, only laminar flows are 
considered in order to remove any uncertainties involv- 
ing the turbulence modeling. 


with 


The equations are written in dimensionless form 
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The subscript re / denotes reference quantities, and for 
convenience the tildes are dropped from the equations. 


Pseudocompressible Formulation 


The system of equations given above is not as 
readily solvable by a fast alternating-direction-implicit. 
(ADI) scheme as the unsteady compressible Navier- 
Stokes equations which are hyperbolic. This difficulty 
can be alleviated by modifying the continuity equation 
as follows (see refs. 11,12): 


1 dp dm 

P dt dxi 


( 1 *) 


The parameter 1//3 is the pseudocompressibility. As 
the solution converges to a steady state, the pseudo- 
compressibility effect approaches zero yielding the in- 
compressible form of the equations. 


For time- dependent problems, a pressure itera- 
tion term is added to the right-hand side of equation 
(If) as follows: 


1 dp dui 1 dp * 

P~dt + Jx 7 ~~ P~df 


( 1 *) 


Here, p* is the value of p at the previous iteration. 
The pressure must be iterated until the divergence-free 
velocity field is obtained within the desired numerical 
accuracy. 


Physical Meaning of the Pseudocompressibility 


The system of equations governing the motion 
of an incompressible flow is elliptic. In the elliptic for- 
mulation, pressure waves propagate with infinite speed. 
However, the system of modified equations given by 
equations (lb) and (If) is hyperbolic. Thus, in the 
present formulation, waves of finite speed are intro- 
duced. The magnitude of the wave speed depends on 
P. In viscous flows, the behavior of the boundary layer 
Is very sensitive to the local presure gradient, espe- 
cially when the boundary layer is separated. For ex- 
ample, if the local pressure gradient oscillates owing 
to the traveling pressure waves of finite speed that are 
generated numerically, the location of the flow separa- 
tion may fluctuate accordingly. This unsteady be- 
havior of the separation will feed back to the pressure 
field, making it difficult to obtain steady-state solu- 
tions. Therefore, the following aspects are of prime 
importance for the success of the present method: 
(1) how the introduced waves of finite speed must 
propagate relative to the viscous effect to simulate 
the incompressible flows, (2) how to choose the value 
of p to achieve computational efficiency, and (3) 
how to guarantee the incompressibility and the ac- 
curacy of the solution. These aspects are analyzed in 
previous papers [11,12], and numerical experiments are 
presented here to demonstrate the characteristics of 
the present procedure. 


Lower Bound of p 

Equation (If) shows that dui/dxi approaches 
zero when P >> 1 . However, since the magnitude 
of p controls the speed of the pressure wave, it plays a 
very important role in determining convergence speed, 
accuracy, and stability. By simplifying the govern- 
ing equation to a one-dimensional linearized form, 
propagation characteristics of the pressure wave can be 
studied [12]. If the velocity in the primary flow direc- 
tion is represented by u, the pseudosound speed c for 
this system of equations is 

c = vW/» (2) 


The dimensionless time scale, tj_, (i.e., the time it takes 
for the upstream propagating wave to travel a distance 
L) becomes 


tl 


L/x„j 
c — u 


(3a) 


The spreading of the viscous effect can be represented 
by the boundary-layer growth rate. Therefore, for 
laminar flows, the time required for the viscous effect 
to spread over a thickness 8 is approximately given by 



(3b) 
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To recover the incompressible phenomena, the physics 
requires that the pressure wave propagates much faster 
than the spreading of vorticity, namely, 

t l <<Te (4) 

From this, a criterion for the lower bound of ,9 is ob- 
tained as follows: 

(5 *» 

When this requirement is not satisfied, the numerical 
procedure can diverge or the rate of convergence be- 
comes so slow that it is practically impossible to ob- 
tain a converged solution. A similar analysis leads to 
a criterion for turbulent flows as follows: 

e»' l +U ! rr)‘&y- 1 (5b > 

Equations (5a) and (5b) set the lower bound of f). The 
upper bound will be discussed in a later section. 

One of the other aspects to be considered is the 
total number of iterations required for the wave to 
travel downstream and back upstream, which is a total 
distance of 2 L. Recalling that the downstream and 
upstream traveling waves propagate with speed (c+ u) 
and (c — u) respectively, the total time required for one 
round trip, zq, can be written as 


T 1 


—J— 2 L 


(6a) 


If a constant time-step is used throughout the com- 
putation, the total number of iterations N% for one 
round trip is 



(6b) 


where A r is the increment of dimensionless time. For 


convergence, enough computing time r has to elapse to 
permit the wave to make at least one round trip. This 
aspect will be demonstrated later using test problems. 


IH. Numerical Algorithm 
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Here, the contravariant velocities, Ui without metric 
normalization, are defined as 


Ui = Lo "I* Liu Lav -j- Law 


( 10 ) 


where 

Lo 1=3 (£•)•> L\ = (£;)n La — (£<)»> La = (£«•)* 

or f for i=l,2, or 3, respectively 

Q 

if = finite difference form of — ,etc. 
h — At — time — step 


Difference Equations 

To accommodate fully three-dimensional geometries, 
the following generalized independent variables are in- 
troduced which transform the physical coordinates into 
general curvilinear coordinates: 

t = t 

t = ((x,jr,z,t) 

, , (7) 

t) = r)(x,y,z,t) 

f = f(*,J/,M) 

Applying this transformation to the governing equa- 
tions (lb) and (If), and combining trapezoidal-rule 
time-differencing and the difference form of the trans- 
formed governing equations, the following equation in 
delta-form is obtained: 


The superscript n denotes the n t h time-step, and the 
viscous terms are given as 

r i - j(v(iV( t i n ^~) (ii) 


Approximate Factorization 

The full viscous terms in Eq. (11) produce non- 
tridiagonal elements in the left-hand side of equa- 
tion (8). Therefore, to implement an approximate- 
factorization scheme, only orthogonal terms are kept 
on the left-hand side. For steady-state solutions, this 
can be done since the left-hand side approaches zero as 
t steady state is approached. For a time-accurate solu- 
tion, this approximation procedure needs to be further 
investigated when a nonorthogonal grid is used. For 
the right-hand side, the full viscous terms may be in- 
cluded. Presently, a nearly orthogonal grid is used, 
and the viscous terms are further simplified to 

A = ^(V^-Vf J/ * = 7 i i m ±- (12 ) 
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After adding smoothing terms to stabilize the com- 
putation, the approximate-factored form of the govern- 
ing equation becomes 

| 7 + \ r +1 6s{Ai n - 'nImS e ) + £<V e A f l- 

[7 + \ J n+1 S v (A2 n - iuI m S„) + e<V, A,]- 

[7 + ^7" +l 5 t (A 3 n - 73/m^) + uV s A ( \(Q n+1 - Q n ) 

= [RHS (7)] - £e[(V,A«) 2 + (V„A n ) 2 + (V f A t ) 2 ]Q" 

(13) 

where and e e are implicit and explicit smoothing 
terms, and 

h — At — time — step 
= Qj — Qi—u — Qk+i Qk 

^Q = (Qi+i-Qi-i)/(2Af) 

8 1,7 8 V Q — [(7fc+l + lk){Qk+i — Qk) 

- (7* + 7k-i)(Qk - Qk-i) 1/[2(A/?) 2 ] 
Analogous terms in the r\ and f directions are defined 
similarly. 

Upper Bound of §_ 

The criteria given by equations (5a) and (5b) fix 
the lower bound of /? for obtaining converged numeri- 
cal solutions. Howevereq, a large value of P makes the 
added pressure term in the continuity equation very 
small, which in turn causes the pseudo-compressible 
formulation to become very stiff [8]. Therefore, the 
upper bound of ft depends on the particular numerical 
algorithm chosen. For example, if an explicit scheme is 
used, the stability condition associated with the pseu- 
dosonnd speed must be satisfied just, as in the case of 
compressible flows. 

The approximate-factorization algorithm used in 
the present study has a different feature to be con- 
sidered. In factoring the finite-difference form of the 
governing equations to obtain equation (13), higher- 
order cross-differencing terms are added to the left- 
hand side of equation (8). The added terms consist of 
the product of A as follows: 

^(J n+1 f5(Ai n 8 v A2 n AQ, 

\ (14a) 

^-(J n+1 ) a S ( A 1 n 8 v A 2 n S { A 3 "AQ, ...etc. 

8 

where 

A Q = Q n+l - Q n (14b) 

These added terms contaminate the momentum equa- 
tions as well as the continuity equation, and there- 
fore must be kept smaller than the original terms 
everywhere in the computational domain; that is, 

o|(^7) 2 M^„A 2 J < (15) 

From the expression given by equation (9), this re- 
quirement leads to the following criterion for t he upper 
bound of p: 

ph < 0(1) (16) 

If A is chosen to be larger than the above range, 
the accuracy of the numerical solution will deteriorate 
significantly and can even induce numerical instability. 


IV. Computed Results — 

Numerical experiments are performed to verify 
both the present algorithm and the criteria for p which 
were based on a one-dimensional linear analysis. The 
magnitude of P is determined by the geometric and 
physical parameters such as the characteristic length 
of the body in the flow direction, Reynolds number, 
and the time-step size used for iteration. Also the 
influence of the p term on convergence, rate and ac- 
curacy is investigated. As a measure of convergence, 
root-mean- square values of A Q (denoted by RMSDQ) 
in equation (13) are monitored. The accuracy of the 
solution can be tested by checking how accurately the 
divergence-free velocity field is attained. The first term 
of A Q in equation (13) can be regarded as a parameter 
showing how well the modified continuity equation (If) 
is satisfied. The root-mean-square value of this quan- 
tity is denoted by RMSCO, which contains smoothing 
terms of the present algorithm. The true incompres- 
sibility is tested by using the root-mean-square value 
of (div u). In the following test problems, these quan- 
tities, namely, RMSDQ, RMSCO, and RMS(div u), 
will be shown to be functions of other parameters. 

Flow through a Channel 

The channel flow is perhaps the simplest internal 
flow test problem, where the pressure wave propagates 
between the inflow and outflow boundaries while the 
viscous effect spreads inward from two walls. The 
coordinate system used for this problem is shown in 
figure la, where velocity vectors for a converged solu- 
tion are also shown. To obtain fully developed velocity 
profiles within a reasonable channel length, a par- 
tially developed boundary-layer profile is imposed at 
the inflow boundary which leaves a flat core region in 
the middle. In figure lb, a converged velocity profile 
is compared with a fully developed velocity profile. 
In this experiment, the reference length and reference 
velocity are 

Xrtj = 2 h = channel height 

Uref “ U a vg 

Tests are performed using the following three cases: 
Casel.L = 20, Re = 1000, At = 0.1 : 0.75 < P < 10 

Casei. L — 30, Re = 1000, At = 0.1 : 1.19 < 0 < 10 

CaseZ. L = 40,Re= 1000, A r = 0.1 : 1.69 < p < 10 

The Reynolds number and the time-step size are kept 

unchanged while the dimensionless channel length is 
varied from 20 to 40. This effectively changes the ratio 
of the time scales required for the pressure waves and 
the vorticity to map the entire flow field. For example, 
the time required for the pressure wave to make one 
round trip from the inlet to the exit boundary will be 
doubled when the channel length is increased from 20 
to 40 while the vorticity spreads to the middle of the 
channel in about the same time. This requires a higher 
wave traveling speed in the algorithm. In the test 
cases, the minimum and the maximum value of p are 
estimated using equations (5a) and (16), respectively. 
In table 1, the number of iterations for one round-trip 
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by the pressure wave is tabulated for various values 
of § which include values outside the recommended 
range. 

In figure 2, RMSDQ histories are shown for the 
above three cases. For all cases, p = 0.1 is much 
smaller than the minimum p required. For L=20, the 
iteration procedure becomes unstable after 340 steps, 
and for L==30 and L=40, the solution seems to be 
converging. 

In figure 3, RMSDQ histories are plotted in terms 
of the number of round trips made by the pressure 
waves, N/N\, between the inlet and exit boundary of 
the channel. In each case, runs with different p are 
plotted at the same interval of the iteration count; that 
is, in figure 3a, the number of iterations using p = 
1 is approximately 3 times greater than the number 
using p = 5 for the pressure wave to reach a same 
distance. These plots show that most variations in 
RMSDQ take place before N / N t reaches 1. By the 
time waves complete two round trips, that is N/Ni = 
2, RMSDQ values are fairly well converged. Using ,8 = 
0.1, the number of iterations to reach N/Ni = 1 is 
4196 for L — 20 and becomes 8391 for L = 40. These 
requirements are not acceptable for practical applicar 
tions even if the solutions converge. Moreover, the ac- 
curacy deteriorates, as will be shown next. Therefore, 
the computation is not continued to reach N/Ni = 2 
for case 3. 

In figure 4, RMS(div u) values are plotted to 
check the accuracy of the converged solutions. When 
the value of p is out of the range specified for the test 
cases, the accuracy of the solution measured by the 
incompressibility deviates from that expected by the 
truncation error analysis of the differencing scheme. 
For example, when p is increased to 50, the incompres- 
sibility is not very well satisfied. Using p = 0.1, the 
accuracy becomes very poor. Since the Mach number 
based on the pseudosound speed in equation (2) be- 
comes very large in this case, the pseudocompressible 
effect becomes significant. This causes the vorticity 
and the pressure wave to interact improperly, failing 
to produce an incompressible effect. This aspect will 
be illustrated further in figure 6. 

In figure 5, RMS(div u) values are plotted in 
terms of N/Ni. These plots also show that about two 
round trips of the pressure waves are required to obtain 
converged solutions. These experiments show that a 
lower value of p within the required range reduces 
the approximate-factorization error at the expense of 
increased iteration. 

In figure 6, RMSCO and RMS(div u) are com- 
pared for case 2. For p — 0.1, which is much smaller 
than required, RMSCO is converging which indicates 
that the modified continuity equation (If), including 
smoothing terms, is reasonably well satisfied. However, 
RMS(div u) deviates from RMSCO considerably be- 
cause of the pseudocompressible effect of the present 


formulation. For p = 1, which is slightly smaller than 
the lower bound, RMS(div u) is converging slowly, even 
though its magnitude is still substantially different 
from that of RMSCO. For p = 2, which is within 
the recommended range, RMSCO and RMS(div u) are 
both converging satisfactorily. 

In figure T, the RMSCO history for all cases is 

shown. 

Flow through a Turn-Around Duct 

To test internal flows further, an annular duct 

O 

with a 180 bend is chosen. This configuration is 
similar to the turn-around duct of the hot-gas manifold 
in the Space Shuttle main engine (SSME). In figure 8, 
the geometry and a laminar solution at Re=l,000 are 
shown, which reveals the formation of a large separated 

O 

bubble after the 180 bend. For this geometry, the 
streamwise length normalized by the duct width is 20. 
Using equation (5a) and (16), the values of p are ob- 
tained for Re = 1,000 and 100 with At = 0.1 

0.17 <<£<10. for Re =1,000 
2.24 « p < 10. for Re = 100 

The test problems presented here were treated using a 
51 x 17 x 21 mesh for half-duct formulation. 

In figures 9a and 9b, the influence of p on the con- 
vergence and the accuracy is shown for Re = 1,000. 
As a measure of incompressibility and the convergence, 
the RMSCO values are shown in figure 9a. As the 
solution converges to the steady state, this quantity 
must approach zero to satisfy the continuity equation. 
Another important quantity to observe is the maxi- 
mum total pressure at any point in the computational 
domain. Here, the total pressure is defined by the fol- 
lowing expression: 

^ Vo Pol 

(-spO 

Pol 

where 

V« — V + ^ + w 2 ) 

and Poi is the total pressure at the inlet. The depen- 
dence of the maximum total pressure on the magnitude 
of p is shown in figure 9b. This experiment clearly 
shows that if p is chosen either too large or too small, 
the solution is not converging. Moreover, if the value 
of p is too large, the approximate-factorization scheme 
is contaminated to the point that the incompressibility 
is not satisfied. 

A similar experiment is done for the case with 
Re — 100, and the results are shown in figures 10a 

and 10b. 

External-Flow Test Problem 

To represent an external flow, the flow past a 
circular cylinder at a Re — 40 was chosen. For ex- 
ternal flows, in which the computational domain ex- 
tends a large distance from the body, the pressure 
waves originating from the body surface propagate into 
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the far-field. Therefore, to obtain the near-field solu- 
tion only-, the distance traveled by the waves and the 
spreading of the vorticity can be approximately the 
same in magnitude. In the present case, if the vis- 
cous region is taken to be approximately two diameters 
away from the body, equation (5a) leads to the follow- 
ing range for /? using A r = 0.1 : 

0.1 << £ < 10 

This indicates that the magnitude of /J is less restrictive 
for external flows in general. 

In figures 11a and lib, the stream-function con- 
tours and the pressure coefficient on the surface are 
shown for a steady-state solution. This solution agrees 
very well with that of Mehta who used a stream 
function and vorticity formulation in two dimensions 
(personal communication from U. B. Mehta, 1983). In 
figure 12, where the history of the pressure drag is 
shown for an impulsively started circular cylinder at 
Re = 40, four different values of /? were compared 
with the time-accurate solution of Mehta. In all cases, 
the values of ft are selected within the suggested range 
above, and the solutions converge quite rapidly. 

The computing time required for the present ex- 
periments was 1.1 x 10 — 4 sec per mesh point per time- 
step on the Cray X-MP computer at NASA Ames 
Research Center. 

V. Concluding Remarks 

This paper presents salient features of the com- 
putational method developed for a three-dimensional, 
incompressible, Navier-Stokes code (INS3D) [11]. 
INS3D has been applied to various geometrically com- 
plex flows, including a major application in analyzing 
the complex flow field in the Space Shuttle main en- 
gine power head. The present study clarifies physical 
and numerical characteristics of the pseudocompres- 
sible approach in simulating incompressible flows. The 
present algorithm has been shown to be very robust 
and accurate if the selection of /? is made according to 
the guidelines presented here. For time-accurate solu- 
tions and also for improvement in the convergence rate, 
an iteration algorithm utilizing equation (lg) is being 
investigated. The effectiveness of this method will be 
published in a future report. 
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Table 1 : Number of iterations required for one round-trip by 
pressure waves between in- and out-flow boundary of 
a channel: Re = 1000 and At = 0.1 
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Figure 8.- Flow through a turn-around duct, (a) Three-dimensional grid, (b) Typical flow pattern with separation (Re = 1 ,000). 
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Figure 9.— Influence of 0 on incompressibility and convergence for 
turn-around duct at Re = 1 ,000. 



NUMBER OF TIME-STEPS 

Figure 10.- Influnce of (3 on incompressibility and convergence for 
turn-around duct at Re = 1,000. 
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Figure 1 1 Steady-state solution for flow over a circular cylinder at 
Re = 40. (a) Stream-function contours, (b) Pressure coefficient on the 
cylinder. 
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